knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
library(tidyverse)
library(ggplot2)
library(ggthemes)
library(kableExtra)
library(mgcv)
library(ggpmisc)
library(mgcViz)
library(viridis)
source('scripts/TSS_mclapply_ie_stn_month.R')
load.Rdata <- function( filename, objname )
{
d1 <- load( filename )
eval( parse( text=paste( objname, "<<- ", d1 ) ) )
}
dat<- as.data.frame(read.delim("data/cal_vert_data_std_light_bloom_for_GAMs.txt"))
month_assign_cfin=data.frame(Month=1:12, TAXON="Cfin") %>% mutate(season=ifelse(Month >6| Month <3, "D", "A"))
month_assign_chyp=data.frame(Month=1:12, TAXON="Chyp") %>% mutate(season=ifelse(Month >3, "D", "A"))
month_assign_cglac=data.frame(Month=1:12, TAXON="Cglac") %>% mutate(season=ifelse(Month %in% c(4,5), "A", "D"))
month_assign<- rbind(month_assign_cfin,month_assign_chyp,month_assign_cglac)
dat<- left_join(dat %>% mutate(fMonth=as.factor(Month),
cum_n_m2 = ind.m2_total * pcum), month_assign) # for simulation
theme_cl<- theme_few()+theme(axis.title.x = element_text(size=14, colour="black"),
axis.text.x = element_text(size=13, colour='black'),
axis.title.y = element_text(size=14, colour="black"),
axis.text.y = element_text(size=13, colour='black'),
legend.title=element_blank(),
panel.border = element_rect(colour = "black"),
axis.ticks = element_line(colour="black"),
strip.text.x = element_text(size=14, face="bold"),
strip.text.y = element_text(size=14, face="bold"))
reducing eps in gam formula with betar distribution has a positive effect on the normality of random effects.
criteria for correction: 1) % of depth sampled < 95% (tow_depth/sta_depth) 2) sta_depth-tow_depth > 10. if sta_depth-tow_depth < 10, then considered as sampled properly.
#to be replaced by Ben package?
library(tidyxl)
ecomon<- readxl::read_excel(sheet="Data","C:/LEHOUX/Projects/Données_copépodes/ACCASP2020/Calfin_Thru_12_30_2019_10m2.xlsx")%>% mutate(perc_sampled= tow_depth/sta_depth *100,
to_correct=ifelse(perc_sampled < 95 & (sta_depth-tow_depth) >10, "Y", "N"),
Zcat = as.factor(ifelse(sta_depth < 200, "0-200",
ifelse(sta_depth >=200 & sta_depth < 300, "200-300",
ifelse(sta_depth >=300 & sta_depth < 500, "300-500",
ifelse(sta_depth>=500 & sta_depth < 1000, "500-1000",
ifelse(sta_depth>=1000, "1000+" , NA)))))))
ecomon %>% group_by(to_correct,Zcat) %>% tally()
## # A tibble: 8 × 3
## # Groups: to_correct [2]
## to_correct Zcat n
## <chr> <fct> <int>
## 1 N 0-200 20502
## 2 N 200-300 627
## 3 N 300-500 3
## 4 Y 0-200 1403
## 5 Y 1000+ 5
## 6 Y 200-300 1563
## 7 Y 300-500 269
## 8 Y 500-1000 15
load("gamms/Cfin_CVIgamms3_4_remove01.RData"); dattax_month <- dattax %>% mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CVIgammsseason_remove01.RData"); dattax_season <- dattax %>% mutate(fMonth=as.factor(Month))
AIC(mod3,mod4, modseason)
## df AIC
## mod3 73.65314 -980.0647
## mod4 47.62006 -2209.9254
## modseason 60.45941 -950.3503
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)
# kbl(tss_table) %>%
# kable_paper(full_width = F)
### Month as continuous
### Season
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod3),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod3, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod4),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_season <- dattax_season %>% mutate(residuals=residuals(modseason),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
dattax_season$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")
print(plot(getViz(mod3), all.terms=T), pages=1)
new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
Zstation = seq(50,500,10),
percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
as.data.frame(predict(mod3, newdata = new_data,
se.fit = TRUE)))
ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
colour = fMonth)) +
geom_line() + theme_cl+
facet_wrap(~ fMonth)
#re<- as.data.frame(predict(mod3, type="terms")) %>% select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl
print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))
print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
print(plot(getViz(modseason), all.terms=T), pages=1)
datsp <- dat %>% filter(TAXON_STADE=="Cfin_CVI")
dat200 <- datsp %>% dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>% dplyr::slice(which.max(ZMax))
dat_zstation <- datsp %>% dplyr::group_by(ID) %>% slice(which.max(ZMax))
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
load("gamms/Cfin_CVgamms3_4_remove01.RData"); dattax_month <- dattax %>% mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CVgammsseason_remove01.RData"); dattax_season <- dattax %>% mutate(fMonth=as.factor(Month))
AIC(mod3,mod4, modseason)
## df AIC
## mod3 237.2257 -4282.196
## mod4 149.0354 -11564.548
## modseason 157.3500 -4231.050
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)
# kbl(tss_table) %>%
# kable_paper(full_width = F)
### Month as continuous
### Season
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod3),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod3, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod4),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_season <- dattax_season %>% mutate(residuals=residuals(modseason),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
dattax_season$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")
print(plot(getViz(mod3), all.terms=T), pages=1)
new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
Zstation = seq(50,500,10),
percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
as.data.frame(predict(mod3, newdata = new_data,
se.fit = TRUE)))
ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
colour = fMonth)) +
geom_line() + theme_cl+
facet_wrap(~ fMonth)
#re<- as.data.frame(predict(mod3, type="terms")) %>% select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl
print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))
season
print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
print(plot(getViz(modseason), all.terms=T), pages=1)
datsp <- dat %>% filter(TAXON_STADE=="Cfin_CV")
dat200 <- datsp %>% dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>% dplyr::slice(which.max(ZMax))
dat_zstation <- datsp %>% dplyr::group_by(ID) %>% slice(which.max(ZMax))
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
load("gamms/Cfin_CIVgamms3_4_remove01.RData"); dattax_month <- dattax %>% mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CIVgammsseason_remove01.RData"); dattax_season <- dattax %>% mutate(fMonth=as.factor(Month))
AIC(mod3,mod4, modseason)
## df AIC
## mod3 185.3930 -3763.047
## mod4 118.6907 -10792.087
## modseason 126.7773 -3780.707
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)
# kbl(tss_table) %>%
# kable_paper(full_width = F)
### Month as continuous
### Season
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod3),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod3, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod4),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_season <- dattax_season %>% mutate(residuals=residuals(modseason),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
dattax_season$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")
print(plot(getViz(mod3), all.terms=T), pages=1)
new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
Zstation = seq(50,500,10),
percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
as.data.frame(predict(mod3, newdata = new_data,
se.fit = TRUE)))
ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
colour = fMonth)) +
geom_line() + theme_cl+
facet_wrap(~ fMonth)
#re<- as.data.frame(predict(mod3, type="terms")) %>% select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl
print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))
season
print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
print(plot(getViz(modseason), all.terms=T), pages=1)
datsp <- dat %>% filter(TAXON_STADE=="Cfin_CIV")
dat200 <- datsp %>% dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>% dplyr::slice(which.max(ZMax))
dat_zstation <- datsp %>% dplyr::group_by(ID) %>% slice(which.max(ZMax))
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
load("gamms/Cfin_CIV_CVIgamms3_4_remove01.RData"); dattax_month <- dattax %>% mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CIV_CVIgammsseason_remove01.RData"); dattax_season <- dattax %>% mutate(fMonth=as.factor(Month))
AIC(mod3,mod4, modseason)
## df AIC
## mod3 276.4344 -4488.609
## mod4 153.7038 -10592.331
## modseason 194.2298 -4322.528
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)
# kbl(tss_table) %>%
# kable_paper(full_width = F)
tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)
kbl(tss_table) %>%
kable_paper(full_width = F)
| ObsPrev_abund | Tresh_opt_e_i | TSS_e_i | TSSsd_e_i | AUC_e_i | corr_r_e_i | p.corr_e_i | intercept_e_i | slope_e_i | r2_e_i | r2sd_e_i |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.65 | NA/NA | NA/NA | NA/NA | NA/NA | 0.83/0.89 | 0/0 | 0.21/0.16 | 0.7/0.73 | 0.65/0.76 | 0.05/0 |
| 0.65 | NA/NA | NA/NA | NA/NA | NA/NA | 0.87/0.86 | 0/0 | 0.26/0.2 | 0.66/0.67 | 0.74/0.7 | 0.07/0 |
| 0.65 | NA/NA | NA/NA | NA/NA | NA/NA | 0.84/0.86 | 0/0 | 0.25/0.21 | 0.63/0.64 | 0.64/0.69 | 0.03/0 |
### Month as continuous
### Season
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod3),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod3, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod4),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_season <- dattax_season %>% mutate(residuals=residuals(modseason),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
dattax_season$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")
print(plot(getViz(mod3), all.terms=T), pages=1)
new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
Zstation = seq(50,500,10),
percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
as.data.frame(predict(mod3, newdata = new_data,
se.fit = TRUE)))
ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
colour = fMonth)) +
geom_line() + theme_cl+
facet_wrap(~ fMonth)
#re<- as.data.frame(predict(mod3, type="terms")) %>% select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl
print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))
season
print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
print(plot(getViz(modseason), all.terms=T), pages=1)
datsp <- dat %>% filter(TAXON_STADE=="Cfin_CIV_CVI")
dat200 <- datsp %>% dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>% dplyr::slice(which.max(ZMax))
dat_zstation <- datsp %>% dplyr::group_by(ID) %>% slice(which.max(ZMax))
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
load("gamms/Cfin_CV_CVIgamms3_4_remove01.RData"); dattax_month <- dattax %>% mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CV_CVIgammsseason_remove01.RData"); dattax_season <- dattax %>% mutate(fMonth=as.factor(Month))
AIC(mod3,mod4, modseason)
## df AIC
## mod3 278.2126 -4480.242
## mod4 153.4178 -10597.444
## modseason 194.1278 -4329.660
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)
# kbl(tss_table) %>%
# kable_paper(full_width = F)
tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")
tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)
kbl(tss_table) %>%
kable_paper(full_width = F)
| ObsPrev_abund | Tresh_opt_e_i | TSS_e_i | TSSsd_e_i | AUC_e_i | corr_r_e_i | p.corr_e_i | intercept_e_i | slope_e_i | r2_e_i | r2sd_e_i |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.65 | NA/NA | NA/NA | NA/NA | NA/NA | 0.83/0.89 | 0/0 | 0.2/0.16 | 0.71/0.73 | 0.66/0.76 | 0.05/0 |
| 0.65 | NA/NA | NA/NA | NA/NA | NA/NA | 0.86/0.86 | 0/0 | 0.25/0.2 | 0.66/0.67 | 0.75/0.7 | 0.08/0 |
| 0.65 | NA/NA | NA/NA | NA/NA | NA/NA | 0.85/0.86 | 0/0 | 0.24/0.21 | 0.64/0.64 | 0.65/0.7 | 0.03/0 |
### Month as continuous
### Season
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod3),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod3, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_month <- dattax_month %>% mutate(residuals=residuals(mod4),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point()+
geom_smooth()
dattax_month$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
dattax_season <- dattax_season %>% mutate(residuals=residuals(modseason),
cat_perc= ifelse(percZ_stn < 0.25, "<25",
ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))
ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
geom_point(aes(col=season))+
geom_smooth()
dattax_season$linpred<- as.vector(predict(mod4, type="link"))
ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
geom_point()+
geom_smooth()
ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
geom_point()+
geom_smooth()
print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")
vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")
print(plot(getViz(mod3), all.terms=T), pages=1)
new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
Zstation = seq(50,500,10),
percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
as.data.frame(predict(mod3, newdata = new_data,
se.fit = TRUE)))
ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
colour = fMonth)) +
geom_line() + theme_cl+
facet_wrap(~ fMonth)
#re<- as.data.frame(predict(mod3, type="terms")) %>% select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+
#scale_color_viridis(option="turbo")+theme_cl
#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl
print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))
season
print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")
print(plot(getViz(modseason), all.terms=T), pages=1)
datsp <- dat %>% filter(TAXON_STADE=="Cfin_CV_CVI")
dat200 <- datsp %>% dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>% dplyr::slice(which.max(ZMax))
dat_zstation <- datsp %>% dplyr::group_by(ID) %>% slice(which.max(ZMax))
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")
dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])
dat_corr <- dat_corr %>% mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
ifelse(Zstation >=300 & Zstation < 500, "300-500",
ifelse(Zstation>=500 & Zstation < 1000, "500-1000",
ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))
ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
geom_smooth(method="lm", se=F)+
geom_abline(slope=1, intercept=0)+
theme_few()+
stat_poly_eq(formula = y~x,
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
scale_x_continuous(trans="log10", name="Observed intergrated abundance")